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Abstract 

We describe an asymptotic approach to gated ionic models of single-cell cardiac 
excitability. It has a form essentially different from the Tikhonov fast-slow form as- 
sumed in standard asymptotic reductions of excitable systems. This is of interest 
since the standard approaches have been previously found inadequate to describe 
phenomena such as the dissipation of cardiac wave fronts and the shape of action 
potential at repolarization. The proposed asymptotic description overcomes these 
deficiencies by allowing, among other non-Tikhonov features, that a dynamical vari- 
able may change its character from fast to slow within a single solution. The general 
asymptotic approach is best demonstrated on an example which should be both sim- 
ple and generic. The classical model of Purkinje fibers (Noble, 1962) has the simplest 
functional form of all cardiac models but according to the current understanding it 
assigns a physiologically incorrect role to the Na current. This leads us to suggest 
an "Archetypal Model" with the simplicity of the Noble model but with a structure 
more typical to contemporary cardiac models. We demonstrate that the Archety- 
pal Model admits a complete asymptotic solution in quadratures. To validate our 
asymptotic approach, we proceed to consider an exactly solvable "caricature" of the 
Archetypal Model and demonstrate that the asymptotic of its exact solution coin- 
cides with the solutions obtained by substituting the "caricature" right-hand sides 
into the asymptotic solution of the generic Archetypal Model. This is necessary, be- 
cause, unlike in standard asymptotic descriptions, no general results exist which can 
guarantee the proximity of the non-Tikhonov asymptotic solutions to the solutions 
of the corresponding detailed ionic model. 
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1 Introduction 

1.1 Physiological and mathematical motivation 

Mechanical activity of the heart is controlled by electrical excitation of cardiac cells, char- 
acterized by "action potentials" (APs) across their membranes. Abnormalities of cardiac 
rhythm are a major public health hazard [46], and great efforts are directed to the math- 
ematical modelling of APs. Investigations at various levels of membrane, cellular and 
myocardial organisation have lead to the development of a large number of detailed ionic 
cardiac models (for reviews, see e.g. [13, 19, 22, 25, 44]) for various types of cardiac cells. 
These detailed models are realistic in the sense that they demonstrate good agreement 
with experimental data available to their authors at the time of creating the model. There 
even exists an optimistic view that with the help of detailed cardiac computational models 
"it will soon be possible to do in silico experiments that would be impossible, difficult or 
unethical in animals or patients" [13]. 

In reality however, detailed ionic models of cardiac excitation are immensely compli- 
cated. Their analytical solution is impossible while their simulations are expensive espe- 
cially in three dimensions. Thus, numerous attempts have been made to construct simpli- 
fied mathematical models of cardiac action potentials (APs); for examples, see [2, 5, 15- 
17, 20, 24, 33, 42]. Detailed ionic cardiac models were initially constructed as variations of 
the spectacularly successful Hodgkin-Huxley model of nerve excitability [21]. In a similar 
way the most simplified cardiac models are often based on the elegant FitzHugh-Nagumo 
two-variable reduction [18, 30] of the Hodgkin-Huxley model. A typical way to create a 
simplified model is either to make an appropriate functional generalization of the FitzHugh- 
Nagumo system, or truncate a detailed ionic models. Then the modellers proceed to fit 
the parameters of their simplified models to reproduce the AP shape or other selected 
properties of the chosen detailed ionic model. 

This modelling approach has an obvious fiaw, as such simplified models are not reliable 
outside the range of phenomena which they have been fitted to reproduce. So it would be 
desirable to derive a simplified model from a detailed model, based on a well defined set 
of verifiable assumptions. One possible way to do that is via asymptotic methods, which 
would utilize small parameters available in the detailed model. 

More importantly, reducing the number of equations in the model often delivers only 
minor reduction of its computational complexity. Paradoxically, the better a simplified 
model is, the better it reproduces another difficult feature of realistic models, their stiff- 
ness. That is, detailed ionic models typically have small parameters which considerably 
complicate their simulation. However it would be natural to try and eliminate those pa- 
rameters by asymptotic methods, and resulting problems without small parameters should 
be much easier for computational study. 

Yet another paradoxical property of ionic models is that increasing the number of phys- 
iological details does not necessarily make the models more reliable, as experimental data 
are not always sufficient for unequivocal identification of model parameters. As argued by 
Cherry and Fenton [12], detailed ionic models of the same types of cells in the same species. 
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but developed by different authors, may disagree significantly. Thus there is a demand in 
modelling practice for models that would be realistic physiologically but independent on 
experimentally unreliable details, i.e. depend on fewer parameters. Again, one possible 
way to create such a model is to simplify a "too detailed" model by eliminating exceeding 
details by some sort of asymptotic procedure. 

So there is more than one serious reason to look more carefully into the small parameters 
in detailed ionic models. The progress in this direction was hampered for some time by 
implicit assumption, induced by the success of the FitzHugh-Nagumo system, that the 
essential properties of cardiac APs can be captured by dynamical systems of the type 

e— = f{x,y), -^ = g{x,y), a; G M^ y G M', e < 1, (1) 

in which some of the dynamic variables are fast (x) while others are slow (y) and where 
e > is a small parameter. The asymptotic structure of (1) is mathematically very 
convenient. Indeed, the presence of the small parameter e at the derivatives allows to 
"dissect" equations (1), in the leading order in e, into a slow-time (degenerate) subsystem 

= /(x,y), ^ = g{x,y), (2) 



and a fast-time subsystem 



dt 



^ = /(^,y), § = 0, (3) 



with T = et, which are much easier to study. It is assumed that the fast subsystem (3) 
has, for all relevant values of v, no attractors but isolated, asymptotically stable equilib- 
ria. There exist a classical theory stating the general conditions which guarantee that, in 
the limit of e — > 0, any finite segment of the graph of solution of the full equations (1) 
approaches uniformly to that of the degenerate subsystem (2) which, in turn, consists of 
slow-motion parts and fast-motion parts separated by junction points. The slow-motion 
parts are on the /-dimensional "slow manifold" determined by the first k equations in (2) 
and the fast-motion parts are arcs of trajectories of (3) within leaves of the "fast foliation" 
y = const [29, p. 173]. A central result of the theory of slow-fast systems is the Tikhonov 
theorem [40]. 

A paradigm for applying the theory of fast-slow systems to biological excitability has 
been laid down by Zeeman [45] . His fundamental idea is a generalization of the FitzHugh- 
Nagumo view and states that the resting state is the unique and globally stable equilibrium 
of the full system, but in the fast subsystem, it is only one of three equilibria. There is an- 
other stable equilibria corresponding to the excited state, and the two stable equihbria are 
separated by an unstable equilibrium which thereby represents the threshold of excitation. 
The repolarization from the excited state, which is an equilibrium in the fast subsystem 
but not in the full system, to the resting state which is the true equilibrium, happens via 
the slow system, and its details, i.e. the shape of the action potential, depends on the 
structure of the slow manifold. 
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Thus, we refer to excitable systems with the asymptotic structure of equations (1) as 
FitzHugh-Nagumo type or Tikhonov-Zeeman systems. 

Extension of this ideology to spatially-extended excitable systems produces a very at- 
tractive and promising asymptotic theory, see e.g. [41] for a review. In this theory, descrip- 
tion of excitation waves is decoupled into description of fast motion of their sharp "fronts" 
and "backs", and description of slow parts of APs outside fronts and backs. After appro- 
priate rescaling, both the fast and the slow subsystems do not have the small parameter 
in them, so their numerical simulation could be implemented efficiently. 

However, this attractive theory was never really applied to cardiac ionic models. It ap- 
pears that FitzHugh-Nagumo type systems fail, in principle, to reproduce some qualitative 
features of ionic models. Here are some examples. 

Features of cardiac excitability: 

1. Slow repolarization. Cardiac APs tend to have very fast upstrokes and much 
slower other phases, including repolarization ( "downstroke" ) . In the asymptotic limit 
e — s> 0, a FitzHugh-Nagumo system of the form (1) with / = 1 will have a fast, i.e. 
duration t ~ e, upstroke of the action potential will have also fast, ~ e downstroke. 
In principle, this can be avoided if the slow manifold, defined by f{x,y) = 0, has a 
cusp singularity with respect to foliation y = const, which is theoretically possible 
if / > 2 [45]. However, our attempts to identify such a cusp singularity in cardiac 
equations have not been successful [11, 38] 

2. Slow subthreshold response. An excitable system reacts to a subthreshold stim- 
ulus by an immediate return to the resting state. In the asymptotic limit e — > 0, 
a FitzHugh-Nagumo system of the form (1) will have both subthreshold return and 
super-threshold upstroke very fast, t ~ e. However, subthreshold return in real cells 
and realistic models has speed comparable to the slow stages of the AP, i.e. much 
slower than the upstroke. 

3. Fast accommodation. If the stimulus is applied not instantly but gradually, then 
the threshold of excitation may increase, and if the perturbation is too slow, the 
system may fail to generate AP altogether (phenomenon known as accommodation). 
FitzHugh-Nagumo systems demonstrate accommodation but only if the time scale 
of the stimulus is comparable to the duration of the AP. In real cells and realistic 
models, accommodation is observed for much faster stimuli, comparable more to the 
upstroke duration than to the AP duration. 

4. Variable peak voltage. The maximum of the AP in a single cell does not, in the 
first approximation, depend on the way the AP has been elicited. Moreover, the 
asymptotic theory of [41] predicts that the maximum of the AP in a propagating 
wave will be the same as in a single cell. However, in real cells and in realistic models 
the maximum of AP does depend on the mode of excitation, and in the propagating 
AP it could be significantly different than in a single cell. 
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5. Front dissipation. The fast accommodation of realistic models has an interesting 
consequence for propagating waves. If excitability ahead of the wave is temporarily 
blocked by a transient process, the wave will expire if the block lasts too long. In a 
FitzHugh-Nagumo type system, this will happen if the duration of the block equals 
the duration of the wave; that is, the wave will cease when its "length", understood 
as the spatial size of the excited zone, has decreased to zero [43]. In contrast, in ionic 
models the front looses its sharpness, "dissipates", much sooner than AP finishes, 
and after that fails to propagate even though AP continues in many cells and the 
excitability of the tissue ahead of the wave has fully recovered [7, 8, 10]. 

In this article, we describe a non-Tikhonov asymptotic approach to cardiac excitation. 
Its purpose is to overcome the limitations of simplified models of FitzHugh-Nagumo type 
discussed above. The approach has already been used in some form in [37] to achieve 
numerically accurate prediction of the front propagation velocity (within 16 %) and its 
profile (within 0.7 mV) for a realistic model of human atrial tissue [14]. The asymptotic 
reduction was sufficiently simple to allow the derivation of an analytical condition for 
propagation block in a re-entrant wave which was in an excellent agreement with results 
of direct numerical simulations of the realistic atrial ionic model. This has been achieved 
by considering only the "fast subsystem" of the full model. 

Here we take a further step and present a complete asymptotic description of a simple 
cardiac excitation model, including both fast subsystem and slow subsystem, i.e. describing 
the whole AP rather than the upstroke only. A brief sketch of this description has been 
outlined in [9]. That sketch has left several questions unanswered, and the purpose of the 
present article is to fill in the gaps. The asymptotic reductions we propose are based on 
a well-defined and verifiable set of assumptions and in this sense they are "derived" from 
a detailed ionic model. We do not include arbitrary fitting parameters which limit the 
model to the reproduction of few hand-picked AP properties. All arbitrariness is restricted 
to choice of small parameters in the model, which is the key stage in any asymptotic 
approach. In our approach, the main small parameter occurs in an essentially different way 
from (1). Consequently, the results of the classical theory of slow-fast systems [29] are not 
applicable. In other words, there is no existing rigorous theory which can guarantee that 
the asymptotic solutions are close to the true solutions of the corresponding detailed ionic 
model. To deal with this issue, we formulate a caricature model which exactly duplicates 
the asymptotic structure of a detailed ionic model but allows exact analytical solution. We 
proceed to apply the asymptotic procedure to this caricature and to compare the exact 
and the asymptotic solutions in order to validate our approach. 

To demonstrate our approach, we need to chose an appropriate gated ionic model. 
Such a selected model must satisfy two criteria. Firstly it should be as simple as possible. 
Secondly, it should have the generic structure similar to all or most of contemporary ionic 
models. The first requirement is best satisfied by the classical Noble model of excitability 
of cardiac Purkinje fibers [31]. This model has spawned generations of models of increasing 
accuracy and complexity up to modern models with more than sixty differential equations 
per single cell [23]. The Noble model, however, does not satisfy the second requirement. 
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Indeed, while that model reproduces the AP of Purkinje fibers in detail, it does not correctly 
refiect their physiology. As the model was constructed before sufficient experimental data 
on the ionic currents became available [32], the inward sodium current was given the dual 
role of generating the upstroke and maintaining the plateau. To avoid this peculiarity and 
to ensure that our asymptotic procedure is sufficiently general, we propose an "Archetypal 
Model" (AM) which has the generic structure of modern cardiac models but keeps the 
functional simplicity of the the Noble model and is identical to it in the asymptotic limit. 
The asymptotic procedure is then demonstrated on the example of this Archetypal Model. 

The structure of the paper is as follows. We conclude the Introduction with a short 
discussion of the procedure of parametric embedding which is an important instrument 
in our work. In section 2, we use a set of numerical observations to postulate a system 
of axioms for a non-Tikhonov parametric embedding of the Noble model. In section 3, 
we describe the AM and discuss its similarities with contemporary ionic models and then 
study the asymptotic limits of the AM and obtain analytical solutions in quadratures in 
these limits. In section 4, we formulate the exactly solvable "caricature" simplification 
of the AM. In section 5, we validate our general asymptotic results by a comparison of 
the limits of the exact solution of the caricature model to its solution in the asymptotic 
limit. We conclude by outlining the most essential features of our approach and discussing 
possibilities for its application. 

In all asymptotics, we restrict consideration to the leading order, with the exception of 
Appendix A which is about a first-order correction. 

Some parts of this work (subsections 1.2, 2.1, 2.3, 3.1, 3.2) appeared in a very brief form 
in [9] and are reproduced here in the full form, with additional details, for completeness 
and convenience of reference in the rest of the article; other parts (subsections 2.2, 2.4, 
sections 4 and 5 and appendices A, B and C ) are entirely new. 

1.2 The procedure of parametric embedding 

The nonlinear problems of physical, chemical and biological applications do not normally 
have parameters which are literally approaching zero within their normal range relevant to 
the application. Hence, a typical practical approach to asymptotic reduction is to identify 
a "small constant", say a, in the model and to replace it by a parameter e, so that the 
original problem corresponds to e = a, whereas the asymptotic formulae are obtained in 
the limit e — > 0. A mathematically equivalent modification of this procedure is based on 
the following 

Definition A parametric embedding with parameter e of a function f{x) is any function 
/(x; e) such that f{x, 1) = f{x) for all x G dom(/). A parametric embedding in the context 
of e — > is called asymptotic embedding. An embedding of a dynamical system corresponds 
to an embedding of its generating vector field or map. 

Thus, the "small constant" a is replaced by ea, and then the original problem corre- 
sponds to e = 1 while the asymptotic analysis is performed in the limit e ^ 0. As long 
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as a is a nonzero constant, the limits e — » and ea are mathematically equivalent. The 
purpose of the above definition is uniformity, especially when the small parameters appear 
in more than one place in the equations. From this perspective, the algorithm of obtaining 
an approximation using, say, a partial sum of an asymptotic series is: the power series in e 
is truncated to a selected number of terms, and then e = 1 substituted in the result. When 
applied formally, this may look counter-intuitive, and yet for reasons explained above, this 
is precisely equivalent to what is always done when "small quantities up to a certain order" 
are taken into account in any asymptotic approach. 

So, asymptotic analysis of a mathematical model by necessity implies introduction of 
artificial small parameters, which is equivalent to drawing a curve in functional space, 
/(x;e), with the only condition that at e = 1 this curve passes through the given point 
f{x). There are infinitely many ways to do this, and the question arises, which of these 
embeddings is "the correct" or "the better" one. Unless the resulting asymptotic series for 
the solutions are converging and the error terms can be estimated, there is no obvious an- 
swer to this question within a purely mathematical context. So the choice should be based 
on practical considerations. The embedding should be such as to allow efficient asymptotic 
analysis. A better embedding should also provide a better quantitative approximation for 
a selected class of solutions to the original problem, and/or preserve better their qualitative 
features of interest. For a given embedding, these aspects can be verified by comparing 
solutions at e = 1 with solutions at smaller e. In terms of the more conventional "small 
constant" a approach, the procedure is to verify if a is indeed small enough to be used for 
asymptotics and try to reduce it further and see how the solutions behave. Note this can 
be done numerically, prior to any analytical work. 

In this paper, we are interested in AP solutions. Typically, we will propose an embed- 
ding and assess its quality by comparing the numerical AP solutions at e = 1 and e = 10~^. 
If the embedding is found reasonable, we proceed to study the limit e — analytically. 



2 . 1 Formulation 

The original Noble [31] model may be simplified by an adiabatic elimination of the super- 
fast m-gate. 



2 The Noble model 



dE_ 

'dt 
dh 

It 



U{E) {h^{E)-h), 



g,{E)ml{E)h + g,{E)n' + gs{E) 



(4b) 





(4c) 
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where 

gi{E) = CAr'gNa {ENa " E) , g2{E) = Cm-^Qk {Ek - E) , 
g^{E) = Cm-' [gNa, {E^a - E) + gK^E) {Ek - E)] , 
gK, (E) = 1.2 exp {{-E - 90)/50) + 0.015 exp {{E + 90) /60) , 
VooiE) = ay{E)/ {ay{E) + /3y{E)) , y = h,n,m, 
fy{E)=ay{E)+(3y{E), y = h,n, 

a (E) = (-^-''^ 3 (E) = (^ + ^) 

'"^^ exp((-E-48)/15)-l' ^""^^ exp((i? + 8)/5)-l' 

a.iE) ^ 0.17 exp ((-E - 90)/20) , ^ ((_^ + 1 ' 

"'^^^^ ^ expTi^i"50VlO)-r - ^(-^ - ^°)/«°) ' 

and 

Cm = 12, (7jva = 400, (7^^ = 1.2, g^a, = 0.14, Ek = -110 and E^va = 40. (6) 

Here E is the transmembrane voltage with values E G [Ek, Ej^^a], and h and n are "gating 
variables" with ranges h E [0,1], G [0,1]; more specifically, h is the inactivation gate 
of the fast sodium current I^a (the first term in the equation for dE/dt) and n is the 
activation gate of the slow potassium current Ik (the second term in the equation for 
dE/dt). Notice that gi{E) > (this represents an "inward current") and (72 < (this 
represents an "outward" current). Our choice for the value of the parameter Ek differs 
from the original value of Ek = —100 used in [31]. This transforms the Noble system [31] 
from a self-oscillatory to an excitable one. Another possibility to achieve this effect is to 
increase the value of the coefficient at the first exponent in gK^, as suggested in [26, 31], 
which leads to similar results [9]. Equations (4) represent a good simplification of the 
Noble model [31] as can be seen by the insignificant difference in the solutions plotted in 
figure 1(a). 



2.2 Naive embedding 

Seeking further simplification, we note that the functions m^{E) and hoo{E) are approx- 
imately stepwise as illustrated in figure 2(a). Thus, as suggested in [6, 8], the crudest 
approximation of (4) is given by 

ml{E)^e{E-Ej, h^{E)^e{Eu-E), (7) 

where E^, E^ satisfy m^{Em) = 1/2 and hoo{Ef/) = 1/2, respectively and 9{-) is the 
Heaviside step function. A similar approximation was done in several earlier simplified 
models, e.g. [16, 17, 20]. They typically took, for simplicity, that E^ = Em- This however 
leads to unsatisfactory description of the front dissipation phenomenon [6]. 
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(a) (b) 




Figure 1: (Color online) (a) Action potential solutions of the original Noble model [31], equations (4), 
and of its "naive approximation", equations (4) with (7). (b) Ionic currents I^a (solid red lines). 
(dashed blue lines) and lout (dotted green lines) in the model (4) (thick lines) and in the model (4) with 
(7) (thin lines). By the tradition accepted in physiology currents that increase voltage are considered 
negative and called "inward" and thus the three currents are defined as iNa = dNa'm^hoo {E — E^a), 
hn = QNai {E - ENa), lout = {gK'n'^ + gKi){E - Ek)- Notice that it is iNa rather than /;„ that is 
the main balance to lout in the full model (4) (thick lines), and there is no such balance for the naively 
simplified model (thin lines). These and all other solutions in the paper are calculated with initial conditions 
£:(0) = -10, h{{)) = 1, n(0) = if not indicated otherwise. 



The naive approximation (7) turns out to be unsuccessful as shown in figure 1(a). The 
AP produced when (7) is used does not have a plateau but returns immediately to the 
resting state. The reason for this behaviour is revealed by an analysis of the individual 
currents which are illustrated in figure 1(b) both for the detailed system (4) and for the 
approximation (7). In the detailed model, the sodium current remains significant during 
the plateau phase, successfully counteracting the potassium current for some time. In the 
approximated model, this current is virtually absent after the initial upstroke. This is due 
to the fact that during the slow decrease of the voltage, both the m and h gates remain 
close to their quasi- stationary values moo, ^oo, and their product W{E) = m^/ioo is exactly 
zero in the approximated model. In contrast, in the detailed model, this product remains 
significant and although it is much smaller than unity, when multiplied by the large factor 
g^a, produces a sodium current which is comparable to the potassium and the leakage 
currents. The current W{E) is called the "window" sodium current, because it runs in 
the region of voltages between Eh and E^ where the gates are supposed to be "almost 
closed" [3]. 

2.3 Axiomatic embedding 

In this section we consider a more elaborate parametric embedding, the asymptotic limit 
of which dissects equations (4) into simpler subsystems. It is based on a number of obser- 
vations of the properties of the Noble model, which will be discussed below, and may be 
formalized in the following 
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Axioms 17 The functions gNa, fh{E), hoo{E) and m^{E) are parametrically embedded 
in the functions gNa{^), fh{E; e), hoo{E; e) and ml^{E; e), e > 0, such that 

1, 



Axiom I. 
Axiom II. 
Axiom III. 

Axiom IV. 

Axiom V. 
Axiom VI. 



9Na[e) -- 
fhiE;e) 



e 9Na, 

e-'fh{E), 



lim ml{E;e) = M{E)e{E-E,^), 

where the function M{E) is close to m^{E) for E > E^ 
lim h^{E;e) = H{E)e{Eh-E) 

where the function H{E) is close to hoo{E) for E < Eh, 
Em. > Eh, 



lim [e-'mliE-e)hooiE;e)j ^ W{E) > 0, 

where W{E) is close to the window current W{E) = m^[E)hoo{E) , 



d 

Axiom VII. \imS{E; e) = lim ( m^(E; e)—h^{E; e 



•dE—' 



0. 



Indeed, the permittivity of the Na current g^a = 400 is large compared to those of the 
other currents gx = 1-2, max((7i^J = 1.8 and g^ai = 0.14 and thus the values of associated 
small constants gx/gNa, Qx = dx, dKi, dNai are of the order of 10^^. This observation is 
formalized in Axiom I by an introduction of a small parameter e which divides gNa- Axiom 
II is postulated on the basis of the observation that the function fh{E) is large compared 
to fn{E) as illustrated in figure 2(b). These functions are reciprocal to the time-scale 
functions of the gates h and n and therefore h is a fast variable while n is slow. The speed 
of h is even comparable to E during the upstroke. These observations are justified in 
[6, 8], where we argued that although in healthy tissue E can be, or at least seems, faster 
than h, there are important applications where the two variables should be considered 
of comparable speed. Axioms III (IV) are suggested by the fact that functions m^(£') 
(/too(-S)) are close to 1 above (below) some switch voltages E^ (Eh), and almost vanish 
otherwise as seen in figure 2(a). The two Axioms do not give a precise definition of the 
functions H{E) and M{E), they only require that these are reasonably close to hoo{E) and 
m^{E) for those values of E where these functions are not small. Here "reasonably close" 
means that replacement of h^{E) with H{E) e{Eh-E) and m^(E) with M{E) e{E-Em) 
should not change significantly the properties of the solution. A possible choice for E^ and 
Eh is so that they satisfy m^^{Em) = 1/2 and hoo{Eh) = 1/2. Axiom V is clearly satisfied 
for equations (4) as shown in figure 2(a). Some simplified models take Em = Eh in similar 
situations [16, 17, 20], however we already know that such simplification affects the "front 
dissipation" feature of realistic models [6]. Axioms III-V have a corollary that 

lmi(rr4(E;e)ME;e)) =0. (8) 
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This indicates that the permittivity of the window current W{E) is small of the order e. 
However, as discussed in subsection 2.2, it is a particular property of the Noble model 
that the window current is finite since it is multiplied by the large factor Qf^fa which is 
of order e~^. To ensure this we postulate Axiom VI. Finally, figure 2(b) demonstrates 
the plausibility of Axiom VII where the graph of the function S{E) = S_(E; 1) is shown. 
f:0020 

Thus, according to Axioms I- VII the parametric embedding of the model (4) is 

e~'giiE) rr^{E- e)h + g^iE) + g:,{E), (9a) 

e-^fh{E) {hJ,E-e)-h), (9b) 

U{E){n^{E)-n). (9c) 

First, we consider this system in the fast time T = t/e. Changing the independent variable 
from t to T, taking the limit e — > and using Axioms III and IV, we obtain the fast-time 

system,-^ 

gi{E)M{E)e{E-E^)h, 

ME) {H{E)e{Eh-h)-h), (10) 
0. 

As intended, the right-hand sides of the first two equations are nonzero, thus we have two 
fast variables, E and h. The slow manifold is the set of equilibria of this system and it is 

^Technically, the dynamic variables iJ, h and n as functions of T ought to be denoted by different letters 
than the same variables as functions oit. We follow the tradition, however, and neglect this subtlety. Later 
we comment on some cases to avoid ambiguity caused by this compromise. 



dE 

dT 
dh 

dt 
dn 

dt 



dE 

dT 
dh 

dT 
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defined by the finite equations 

M{E)e{E - Em)h = 0, (11) 
H{E)e{Eh-E)-h = 0, (12) 

since gi{E) > 0, fh{E) > in the physiological range of the voltage, E G [EK,ENa]- 
Substitution of (12) into (11) with account of Axiom V turns (11) into an identity as 
the product of the two Heaviside functions vanishes. Thus, we have a co dimension- 1 slow 
manifold, defined by equation (12). This is a non-Tikhonov feature since in Tikhonov 
systems [29] the codimension of the slow manifold is equal to the number of fast variables. 
This peculiar feature results from (11) becoming an identity if (12) is satisfied, which, in 
turn, is due to a near-perfect switch behaviour of hoo{E) and m^{E), becoming perfect 
switches in the limit e — > 0. A consequence of this feature is that the equilibria of the fast 
system are not isolated. Therefore, Tikhonov's theorem [40] is not applicable as it requires 
asymptotic stability of equilibria of the fast system which does take place here. The 
fact that our parametric embedding is a non-Tikhonov one was already obvious from the 
dependence of (9) on the small parameter which is not of the form allowed by Tikhonov's 
theorem, namely the large parameter in the first equation multiplies only one of the 
terms in the right-hand side, rather than the whole right-hand side. 

Since Tikhonov's theorem cannot be used to describe the slow motion in the standard 
way, we discuss it in more detail, however without attempts of rigorous treatment. We 
consider system (9) in the original (slow) time, and restrict our attention to trajectories 
near the slow manifold, i.e. ones for which h ~ hoo{E; 0). By re-arranging equation (9b), 
and assuming that when moving along the slow manifold the derivative dh/dt is of order 
unity (this assumption is confirmed by the following result), we obtain 

h = h^{E- e) - = M^; e) + 0(e). (13) 

From here we deduce that h = hoo{E; e) + 0{e). Differentiating this with respect to time 
and substituting the result back into the right-hand side of (13), we obtain 

^-^M-J^^'-t'^^Oi.^ (14) 

which substituted in equation (9a) yields 

^ = e^'giiE) m^{E; e) h^{E; e) - |||| S{E; ^ + 0(6) + g^iE) + g,{E), (15) 

where the function S_{E] e) is defined in Axiom VII. Using Axiom VI for the first term and 
Axiom VII for the second term and taking the limit e — ^ we obtain the slow-time system, 

^ = g^{E) W{E) + g^{E) + g:,{E), 

h = H{E)e{Eh-E), (16) 
— = /„(E) {n^{E)-n). 
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Figure 3: (Color online) Solution {E,h,n) of the parametric embedding (9) with explicit expressions 
(17a) and (17b) with p = q = 1, r = 0.5 for e — 1 (thick lines, i.e. equivalent to the authentic model (4)) 
and e = lO^'^ (corresponding thin solid lines) (a) in slow time t G [0, 600], (b) in fast time T G [0, 1]. 



This is a system of two differential equations for the slow variables E and n plus a finite 
equation for h defining the slow manifold. Note that we have two slow variables in agree- 
ment with the two dimensions of the slow manifold. Thus, E is both a fast and a slow 
variable which is yet another non-Tikhonov feature. f:0030 

2.4 Explicit embedding 

To show that Axioms I- VII are consistent and usable we need to demonstrate that there 
exist embedding functions g^a, fh, h^o and which satisfy these axioms. The first two 
functions are already defined by Axioms I and II. Thus, in this section we suggest an 
explicit dependence of and hoc on e which satisfies Axioms III-VI, namely 

mL(E; e) = M{E) (e^ e{Eu - E) + e{E - E,) e{E^ - E) + e{E - EJ) , (17a) 
h^{E; e) = H{E) {e{Eh - E) + e''^' e{E - En) e{E^ - E) + e'^ e{E - Em)) , (17b) 

where p G [l,+oo), q e [l,+oo), r G (0,1) and, of course, we must have M{E) = 
ml^{E), H{E) = hoo{E) to ensure these functions coincide with rn?^{E) and hoo{E) at 
e = 1. It straightforward to verify, that Axioms III-VI are then satisfied. To verify Axiom 
VII is more complicated: obviously, the above derivation of the slow system (16) is not 
technically valid for discontinuous hoo as defined by (17b), and Axiom VII does not make 
sense literally. Still, it will be satisfied if we assume that 6{E — Eh)6{E — Em) = 0, 
i.e. infinity times zero aX E = E^ equals zero. If p = g = 1, we have the asymptotic 
window current exactly equal to the window current of the Noble model, W{E) = W{E); 
if p > 1 and q > 1, then the asymptotic window current is a cut-off version of the original, 
W{E) = W{E) e{E - Eh) e{Em - E). The adequacy of this embedding for p = g = 1, 
r = 0.5 is demonstrated in figure 3. 

The explicit embedding (17) is rather complicated. Formally, there are infinitely many 
embeddings; e.g. equations (17a) and (17b) define a three-parameter family of embeddings. 
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(a) (b) 




Figure 4: (Color online) (a) Solutions {E, h, n) of the Noble model (4) (thick lines) and of the Archetypal 
Model (19) (corresponding thin solid lines), (b) Ionic currents iNa (solid red lines), (dashed blue lines) 
and lout (dotted green lines) in the model of Courtemanche et al. [14] (thick lines) and in the Archetypal 
Model (19) (thin lines). The /ATa-curves of the latter two models overlap. The definition of the currents 
is given in the caption of figure 1. Note that in both models, lout during the most of the AP is mainly 
balanced by /i„ but not iNa- 

all satisfying the Axioms and all leading to the same fast and slow systems and having 
the same asymptotic properties. And it is not possible to infer from the original problem, 
which of the embeddings is "correct" . The discontinuity of hoo and is also a cause for 
concern. It is true that their limits in e ^ have to be discontinuous, or at least non- 
analytical, according to Axioms III and IV, but with (17a) and (17b) they are discontinuous 
already for any e 7^ 1, which is inconvenient even in this heuristic treatment. Moreover, 
this is likely to cause serious technical difficulties in any attempt of rigorous treatment of 
the problem. This difficulty can be overcome, for instance, by using in (17a) and (17b), 
instead of Heaviside functions, functions which are smooth for all e > and only tend to 
Heaviside functions in the limit e = 0, e.g. 

^(x;e) = l/(l + e-^/^). (18) 

Indeed, in that case the peak value of \dhoo{E]e)/dE\ is attained at E = Eh and is 
(4e)~^/ioo(-E'/i) in the leading order, whereas the value of m^{Eh] e) is exponentially small 
in e, thus their product is uniformly small and Axiom VII is satisfied. 

3 The Archetypal Model 

3.1 Formulation and parametric embedding 

f:0040 

The complications arising in the construction of a possible embedding of the Noble 
model (4), as discussed in the preceding section, are not essential. In fact, they arise 
only due to the fact that insufficient experimental data was available at the time of the 
construction of the Noble model and in particular the existence of a Ca current was not yet 
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discovered. Thus, the Na current was given a dual role to produce the upstroke and to keep 
voltage elevated during the long plateau stage leading to the large Na window illustrated 
in figure 1(b). Such a large window current is not present in other cardiac models. This is 
demonstrated in figure 4(b) in the case of the currently-accepted detailed ionic model of 
human atrial cells of Courtemanche et al. [14]. Bearing in mind the possibility of extending 
our present analysis and results to models of other types of cardiac cells, we propose a 
modified version of the Noble model (4). It is similar in many aspects to (4) in but has a 
more "generic cardiac" structure and admits a straightforward asymptotic embedding: 

dE 

— = giiE) M{E) e{E - E„,) h + gi{E) W{E) + g^i^E) + g^{E) 

= gi{E) M{E) e{E - Em) h + g^iE) + G{E), 

dh (19) 

— = U{E){H{E)e{E^-E)-h), 

ci?7/ 

— = /„(E) (noo(^)-n), 

where the functions moo, ^oo, fhi fn, di, 92, 93 are defined as in equations (4), and M[E) = 
ml{E), H{E) = h^E), G{E) = g^{E)W{E) + g^{E), W{E) = ml{E)h^{E). In this 
formulation, the perfect-switch behaviour of the Na gates is represented by the Heaviside 
functions multiplying M and H. The deviation from the perfect-switch behaviour, due to 
the window current component W = m^/ioo, is separated from the fast Na dynamics and 
appears as an additional "time- independent" current with a role similar to S's (-£"). 

We shall call (19) the "Archetypal Model" (AM). The AP of this model is very similar to 
the AP of the Noble model (4) as shown in figure 4(a). Its simplest asymptotic embedding 
consistent with Axioms I-VII can be done with e introduced linearly and only in two places, 

dE 

— = e-'gi{E) M{E) e{E - E^) h + g,{E) W{E) + g^iE) + g;{E), 

^ = e-'UE) {H{E) 9{E, - E) - h) , (20) 

-^ = fn{E) ME)-n). 

The quality of this embedding is illustrated in figure 5. The limit of the fast-time system 
is 

^ = g,(E)MiE)9{E-Em)h, 

^ = h{E) iHiE)eiE,-E)-h), (21) 

^ = 
dT ' 

and the limit of the slow-time system is 



^ = gi{E) W{E) + g^{E) + g:,{E), (22a) 
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Figure 5: (Color online) Solution {E, h,n) of the parametric embedding (20) for e = 1 (thick lines, i.e. 
equivalent to the Archetypal Model (19)) and e ~ 10^^ (corresponding thin solid lines) (a) in slow time 
t e [0, 600], (b) in fast time T G [0, 1]. 

h = H{E)e{Eh-E), (22b) 

^ = fn{E) ME)-n). (22c) 
at 

These systems coincide with (10) and (16), if W{E) = W{E). Their solutions and phase 
portraits are shown in figure 5 and figure 6, respectively. 
f:0050 

The AM (19) produces an AP similar to the Noble model (4). In fact, the agreement 
between the two can be improved further as demonstrated in Appendix A. Moreover, the 
AM and the Noble model have the same asymptotic limits. The AM has two advantages. 
Firstly, the relevant asymptotic embedding is much simpler: the right-hand sides of (20) 
linearly depend on e~^, and it appears only in two places. Secondly and more importantly, 
judging from figure 4(b), the asymptotic properties of the fast Na current in modern 
detailed models are likely to be similar to those of the AM (19) rather than to those of 
the Noble model (4). Therefore, we adopt the AM (19) as the example on which the 
non-Tikhonov asymptotic procedure is demonstrated. 



3.2 Asymptotic analysis 

f:0060 

3.2.1 The fast-time subsystem 

The fast-time subsystem (21) of the AM governs the evolution of the two fast variables E 
and /i on a time scale t ~ e or equivalently T ~ 1. Its solution does not depend on the 
third variable n, which is slow and thus stays close to its initial value no on this time scale. 
System (21) admits exact analytical solution. If Eq < Em the solution of the initial- 
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Figure 6: (Color online) Phase portraits of (a) the fast system (21) and (b) the slow system (22) of the 
Archetypal Model. Blue dashed lines represent vertical isoclines dh/dt = in (a) and dn/dt = in (b). 
Solid red lines and the cross-hatched region in (a) represent horizontal isoclines dE/dt = 0. Thin dotted 
black lines with attched arrows represent trajectories. The thick dotted green line corresponds to the AP 
shown in figure 5(a). The letters A-F designate feature points of the solution. Notice that the blue set in 
(a) is a subset of the red set, so it is a continuous set of equilibria in the fast subsystem. 



value problem with -^(0) = Eq, h{0) = Hq is 

E = Eq, 

h = HiEo) e{En - Eo) + {h - H{Eo) 9{Eh - Eo)) exp (-A(^o) T) . 
If Eq > Em the solution of the same initial-value problem is given in quadratures by 

h = hQ + J{EQ)-J{E), 

E 

drj 



Eo 



gi{v)Miv) {ho + Jm-J{r,)) 



where 



J(r/) 



drj. 



(23) 

(24a) 
(24b) 

(25) 



Note that each trajectory of the fast-time subsystem (21) crosses the slow manifold 
(22b) only once as shown in figure 6(a). Thus, singularity of the slow manifold with respect 
to the trajectories of the fast system, which in Zeeman's [45] interpretation of Tikhonov 
systems determines the threshold of excitation, does not exist in the AM. To shed light on 
the threshold properties of our fast-time subsystem, let us consider the maximal overshoot 
voltage Eoo = E{+oo) in the system as a function of the initial conditions. Taking into 
account that h{+oo) = we can use (24a) to find E^o as a solution of the finite equation 
J (Eoo) = J{Eq) + ho, provided that Eq > E^ and, of course, E^ = Eq, otherwise. Thus, 
the function £'oo(£'o, /iq) has a discontinuity along the line {(£'o,/io)} = {Em} x (0, 1] as 
shown in figure 6(a). This discontinuity is the mathematical equivalent of the physiological 
notion of excitability: the upstroke does take place if and only if Eq is above the threshold. 
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Em- Note that excitability in Tikhonov-Zeeman systems has this property, too, but is 
related to unstable branches of the slow manifold, rather than discontinuities in the fast 
flow. 

The numerical values of parameters of the Noble model from which our Archetypal 
Model descends, are such that fh{E)/gi{E) is a uniformly small function of E in the 
physiological range. That is, condition ho = J(-Eoo) — J{,Eq) ~ 1 requires relatively large 
values of the quadrature (25), the integrand of which is generally a relatively small quantity. 
This is achievable if the integration interval goes close to the pole of the integrand, i.e. 
zero of M{ri). In other words, the noted smallness of fh/gi necessitates that E^o ~ Ej^a 
for typical Eq and ho. A more detailed and formal analysis of this aspect is given in 
Appendix B. 



3.2.2 The slow-time subsystem 

The slow-time subsystem (22) of the AM governs the evolution of the two slow variables 
E and n on a time scale t ~ 1 or T ~ e~^. In fact, equations (22a) and (22c) form a 
closed subsystem. Gate h described by (22b) can be evaluated if E{t) is known but does 
not affect the dynamics of other variables. The slow-time system has been studied by 
Krinsky and Kokoz [27]. The slow-time subsystem (22) is, in turn, a fast-slow system, 
this time in the classical Tikhonov sense. This is illustrated by figure 7(a) which compares 
the characteristic time-scale functions of the dynamical variables E and n, defined as 
Ty = \dy/dy\~^ for y = E,n. It can be seen that £^ is a fast variable and n is a slow 
variable, which motivates the introduction of a second small parameter e2 > in the 
following standard Tikhonov way, 

H E 

^ = G{E)+g^{E)n\ (26a) 
at 

— = e2fn{E) {n^{E)-n). (26b) 



f:0080 

System (26) is a fast-time system. In the limit €2 ^ the lines n = n(to) = const form 
the leaves of the fast foliation. On every such leaf the solution for the voltage is given by 

E(t) ^ 

I G{rt)+gMn{tor ^^^^ 

E{to) 

The corresponding slow-time system is obtained by the change of variable ^2 = £2^ and in 
the limit €2 — > becomes 

= G(E) +^2(^)n^ (28a) 

^ = fn{E) ME)-n). (28b) 
at2 
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Figure 7: (Color online) (a) Characteristic time-scale functions Ty ~ \dy/dy\ ^ for variables y ~ E,n 
of the slow-time system (22) corresponding to the trajectory B-C-D-E-F in figure 5. (b - c) A solution 
{E, n) of the parametric embedding (26) for £2 = 1 (thick lines, i.e. equivalent to the system (22)) and 
£2 = 10^"^ (corresponding thin solid lines) (b) in slow time t € [0,600], (c) in slow time t g [0, 15] and (d) 
in fast time T e [0, 15]. 



Equation (28a) defines the super-slow manifold, 

n = ^^{E) = {-G{E)/g2{E)f\ (29) 

and equation (28b) describes the motion along this manifold. As illustrated in figure 6(b) 
the super-slow manifold is split into two parts by the condition > 0: the "diastolic" 
branch E G {—oo,Ei] and the "systolic" branch for E G [E2,E3], where Ei fii —95.75, 
E2 ~ —61.81 and E^ ^ 1.86 are roots of the equation G{E) = 0. The stability of the fast 
equilibria is determined by the sign of dE/dE which coincides with the sign of N''{E) = 
dM /dE: the stable branches of the super-slow manifold correspond to regions in {n,E) 
plane where its graph has a negative slope, i.e. N''{E) < 0. These are the regions of the 
entire diastolic branch and the upper part of the systolic branch, corresponding to G 
{E^, E's], where E^ k, —17.05 is the root of the equation M^E/) = 0. These considerations 
determine the excitability properties in terms of the super-slow subsystem (28). As seen in 
figure 6(b) the threshold of excitability is E2 since a trajectory starting from -E'(to) > E2 
will be repelled by the lower systolic branch and attracted by the upper one, thus making 
a relatively large excursion. This will be followed by a slow movement along the upper 
systolic branch and a jump to the diastolic branch at E^,. On every monotonic branch of 
the super-slow manifold the equation (28a) can, in principle, be resolved with respect to E 
to give E = ^/'^{n) and the result can be substituted into (28b) leading to the following 
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quadrature solution of equations (28), 



n 



fn{M-\v)) (noo(Ar-l(z/))-Z/)' 



(30) 



where t2,c is the slow time of the beginning of this asymptotic stage. Alternatively, we 
may use, as suggested in [31], E rather than as a coordinate on the super-slow manifold. 
To do that we substitute n = Af{E) into the second equation 



where from solution of the slow-time system in quadratures follows without the need to 
invert the function M{E). In particular, the time between the points (i?, h) = {E^, 0) and 
{E^, A^*), i.e. the duration of the plateau of an AP starting from n = 0, is 



where = ^/{E^) ^ 0.6512 for the standard parameter values. 

This completes a brief overview of the leading-order asymptotics of the Archetypal 
Model. An inventory of resulting formulas describing all the stages of typical solutions 
corresponding to various sorts of initial conditions, is given in Appendix C. 

4 The caricature model 
4.1 Motivation 

As pointed out in the Introduction and throughout the article, the asymptotic structure 
of the embeddings of both the Noble and the Archetypal Model is non-Tikhonov and the 
results of the classical theory of slow-fast systems are not applicable. The fact that we 
are able to demonstrate a good numerical agreement in figures 3, 5 and 7 is reassuring 
but far from reliable since the numerical solutions are, by their nature, only approximate. 
Developing an alternative to the classical slow-fast asymptotic theory is beyond the scope 
ot this paper. Instead, in this section we propose and study a simple caricature model of 
cardiac excitation. One can think of the caricature model as a detailed ionic model which 
allows an exact solution and which has been embedded in a non-Tikhonov system. Once 
exact solutions are available, their the proximity to asymptotic solutions can be proved in 
this particular case. We have to emphasize that, the caricature is not "derived" from the 
Noble or from the Archetypal Model. They are, once again, used only as starting point for 
their simple functional forms. The caricature is constructed so that it has an asymptotic 



dE ^ UE) {n^{E)-M{E)) 
dt2 U'{E) 




(31) 
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structure identical to that of the AM of section 3 and differs from it in that the functions in 
its the right-hand side are chosen so as to make possible to (a) evaluate the quadratures of 
section 3.2 and thus obtain explicit asymptotic solutions and (b) to go a step further and 
find an exact analytical solution of the simple model. We then proceed to to demonstrate 
that the appropriate limits of the exact analytical solution of this caricature model coincide 
with its solution in the asymptotic limits as given by quadratures (23), (24), (27), (30). 
Such an effort is still not a proof of our asymptotic analysis in general, but makes a good 
justification. 

4.2 Formulation 

In order to formulate a caricature of the AM (20) we keep the asymptotic structure of the 
latter and replace ("approximate") the functions forming its right-hand side with simpler 
ones. We make the following simplifications. 

(a) We replace the functions M{E) and H{E) by unity. Thus, the function m^{E) is 
replaced by 6{E — Em) and the function hoo{E) by 9{Eh — E). Analogously, we 
replace the function noo{E) by 9{E — En). The values of the constants Em, Eh and 
En will be discussed in item (e) below. 

(b) We replace the function fh{E) by the constant = 1/2. Analogously, we replace 
the function fn{E) by the constant F„ = 1/270. 

(c) We replace the function G{E) by the continuous piecewise linear function, 

r h{E,~E), Ee(-oo,Et), 
G{E) = l k2{E-E2), Eem,E,), (32) 
[ hiEs-E), E e [E,,+oo), 

where the constants ki = 0.075, k2 = 1/25, k^ = 1/10 while the constants Ei = 
—280/3, E2 = —55 and £^3 = 1 are close to the roots of G{E) (which are ~ —95.75, 
—61.81 and 1.86 respectively; we prefer to avoid multi-digit decimal values). The 
constants E^ = —15 and E^ = —80 are determined by the intersection points of the 
three linear pieces of G{E). Note that at E^ the function G{E) reaches its local 
maximum and thus £"* corresponds to the point on the systolic branch of the super- 
slow manifold which forms the boundary between its stable and unstable parts. 

(d) We replace the function g2{E) by 020 {E — Ej), where g2 = —9. 

(e) Finally, we set Em = E^, E^ = En = E-^. A more accurate approximation would 
be to choose these constants as the solutions of the equations m^^{Em) = 1/2, 
hoo{Eh) = 1/2, and n^{En) = 1/2, respectively. However, this would lead to a 
caricature model, the right-hand side of which would consist of six pieces instead of 
only three as with the present choice. Obviously, this is not a principal difficulty but 
would be rather technically inconvenient. Moreover, the present choice preserves the 
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(a) (b) 




Figure 8: (Color online) (a - b) The caricature approximations (thin lines) of the right-hand side 
functions of the Archetypal Model (19) (thick lines). 



relatively good quantitative agreement with the original models without the addi- 
tional complications. 

A justification for the simplifications of the functions m^(£') and h^{E) can be found in 
figure 2(a) while the rest of the replacements are illustrated in figure 8. 
f:0091 

On these grounds, we postulate the following caricature of the AM (19), 

^ = e-^Gjva (^jva - E) e{E -E,)h + gle{E - E^) + G{E), (33a) 
^ = e-^F, {e{E^ -E)-h), (33b) 

Q77/ 

— = e2F„ {e{E-E^)-n), (33c) 

where we have included the artificial small parameters e and €2 to ensure the correct 
asymptotic structure of the system. 



4.3 Exact solution 

It is possible to obtain an exact analytical solution of the caricature model (33). Indeed, 
equations (33b) and (33c) are separable and simple enough to be easily solvable. After their 
solutions are substituted in equation (33a) it, too, becomes a readily-solvable first-order 
linear ODE. Therefore, assuming the initial conditions, 

E{0) = Eo> E,, h{0) = 1, n{0) = 0, (34) 

of a fast-upstroke AP and natural continuity conditions at the ends of the three intervals 
separated by E-^ and E^ or, equivalently, at the ends of the corresponding time intervals 
t G [0, t^],tE [t*, t-f] and t G [tj, 00), system (33) has the following exact analytical solution, 

[l-exp(-e2F„t), tG[0,tt] 
n{t) = < (35a) 



(exp(e2-F„,t|) - 1) exp(-e2F„t), t G [tf, 00] 
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h{t) 



1 - (1 + exp{Fht^/e)) exp{-Fht/e), te[t^, oo] 

GNa 



(35b) 



exp 



exp 



Fht 

€ 



ht X 



Eo exp 1^ - 

1=0 
GNaENa 



G 



Na 



U 



Fh 
4 
/ 

Fh 
e 



- hEsui-ks^t) 
m((4-0 esF,- A;3,t) 
- k3,t 



(35c) 



E(t) = (E, - w{t,)) exp (A;2 (t - t,)) + M;(t), 
E{t) = (^t - ^i) exp {-ki{t - tt)) + El, 



where 



m(x, t) 



e 



4 



w{t)^E,-9',j2(-^y( 
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Na 



Fh Fh 

A\ e xp (-/ e2F„^) 
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r 



xe G 



Fh Fh 



exp 



t G [0, t,] 
t E [tf, cxd] 



(36a) 
(36b) 



and r(a, x) is the upper incomplete gamma function, T{a, x) = z°'~^e~^ dz for Re(a) > 
and r(a + l,x) = ar{a,x) + x"" as defined in [1]. The exact analytical solution is 
plotted in figure 9 where it is compared with the numerical solutions of the AM (19) and 
the internal boundary points E-^ and E^ are also indicated. The parameters and tj- can 
be found as solutions of 



E{t* 



E^ 



E, 



(37) 



Both equations are transcendental and not solvable analytically. It is possible to solve them 
by using perturbation expansions about the small parameters e and €2. This, however, will 
contribute little beyond its value as a technical exercise and will lead us away from the 
main point of this part, which is to validate the asymptotic solutions for E{t), h{t) and 
n{t). So here we omit these formulae, assuming that t^, and are known where they are 
needed. Numerically, for the standard values of parameters and Eq = —10, we obtain 
^ 292.815 and ^ 345.241. f:0092 
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Figure 9: (Color online) The numerical solution of the embedding of the Archetypal Model (20) (thick 
lines) in comparison with the analytical solution (35) of the caricature model (33) (corresponding thin 
lines) (a) in slow time t G [0, 600], for e = 1, 62 = 1, and (b) in fast time T G [0, 1], for e = 10"^, £2 = 10~^. 



5 Validation of the general asymptotic analysis 

In the following we consider the stages of a normal fast-upstroke AP. For each stage we (a) 
evaluate the appropriate quadratures of section 3.2 to obtain explicit asymptotic solutions 
for this stage and (b) evaluate the limit of the exact analytical solution (35) as the small 
parameters tend to zero as appropriate for the same stage. The asymptotic theory will be 
validated if the results from these two steps are identical. 



5.1 Fast upstroke, stage A— B 

Here and below, letters A-F refer to the labels in figure 5(a). The fast upstroke occurs 
during the interval t G [0,^^] C [0,t*], where is ^ but t^/e — ^ 00 as e — > 0. During 
this stage E and the /i-gate change together on a time scale ~ e from the point {E, h) = 
{Eq, 1) to {Eb, 0) while n is a slow variable and remains approximately at its initial value 
n ^ Uq = 0. 

(a) Asymptotic solution. Asymptotically, this stage is described by quadratures (24,25). 
Substituting there the functions defined in section 4.2, we obtain 

(38a) 
(38b) 

Solving (38b) for the voltage E and substituting the result in (38a), we arrive at an explicit 
asymptotic solution. 




E = E^,~ {Ej,, - Eo) exp (e"^'^^ - l)) , 



(39) 
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The maximal overshoot voltage Eb is obtained as the fast time T tends to infinity, 
Eb = Eoo= lim E(T) = E^a - {E^a - Eq) exp ' 



(40) 



Since exp{—GNa/ Eh) ~ 10~^^, the maximal overshoot voltage Eb is extremely close to 
E^a and h{tB) = h{T oo) = as expected. 

(b) Limit of the exact solution. We change to fast time, T = t/e, and take the limit 
of (35b) and (35c) as e tends to zero, for fixed T. Since linir(a, x) = exp(— x), it follows 

that limM(x, eT) = 0, for x ~ 1 and therefore 



e->0 



limM(— ^3, eT) = 0, 

e— >0 

limM((4-/)e2F„-A;3,eT) = 0. 



Further, 



u 



fe, eT 



G 



Na 



exp 



G 



Na 



— exp 



G 



— e ^ 



i 

Therefore, taking the limit e ^ of expression E in (35c) we obtain,^ 

hmE(eT) = E^a - {E^, - E,) exp (e"^'^^ -I)] 

lim/i(eT) =6"^'^^. 



(41) 



Since equations (39) and (41) are identical, the asymptotic theory of the fast upstroke of 
the AP is validated. 



5.2 Post-overshoot drop of the voltage, stage B— C 

This corresponds to the time interval t G [tB,tc] C [0,t*], where tc — > +oo but e2tc — 
as e2 — > 0. We keep in some of the formulae rather than replacing it with its limit, 
lim^^otB = 0, as a symbolic reminder of the beginning of this asymptotic stage. The 
voltage E decreases on a time scale ~ 1 towards the stable systolic branch, the h-gaie 
remains at /i ^ and the slow n-gate also remains approximately unchanged at n ^ rig = 0. 
(a) Asymptotic solution. The asymptotics of this stage are given by quadrature (27) 
which, with account of the approximations of section 4.2 and of the fact that asymptotically 
n = 0, evaluates to 



^Here we stress that E and h are considered as functions of t rather than T. 
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Solving for E(t), the explicit asymptotic solution for this stage is, 

E{t) = E; + {Eb - E;)e-''^'-''^\ (43) 

(b) Limit of the exact solution. We are still within the interval t G [t_B,tc] C [0,t^] 
1 

and so use E{t) of (35c). Taking the limit of (35a) and (35b) as e and €2 tend to zero, we 
obtain 

lim n{t) = 0, (44a) 

£2— >0 

lim h{t) = 0, (44b) 

in agreement with the asymptotic analysis. 

1 

Now we consider the limit of E[t) given by (35c) for e ^ and €2 — > at a fixed 

t G (0,t^,). In the limit 62 — > the terms u ((4 — /) e2-Fn ~ ^3?^) become independent of 

1 

the index / and the sum over / in the expression for E{t) vanishes since the binomial 
coefficients cancel each other. For the remaining upper incomplete gamma functions, we 
use the recurrence relation, 

r(a + 1, x) = a r(a, x) + x" e'"^ 
and the following asymptotics in the limit a \ 0, 

r(-a,x) = Ei(l,x) + 0(a), 

T{-a,Aexp{-B/a)) = -a'^l-e^)+0{l), ^ ' 

for fixed x, A, B > and where Ei(z/, x) = exp(— X2;) d^;, Re(a;) > is the exponen- 

tial integral. Hence, as e — > 0, we have 



M--fc3,tU7^(expf-%)-ll. (47) 



n(-fc3,t)^fc3~Ml-e'^*)' (46) 

— -h,t \ ^ — 

1 

Substituting these results into the expression for E in (35c) and taking the limits e — >^ 

and €2 — > 0, we get 

lim^ E{t) = + (^Na - {ENa " ^o) exp - - E^ e-'^^*, (48) 

which with account of (40), coincides with the asymptotic expression (43) where is taken 
at its limit, Ib = 0. Thus the asymptotic procedure is validated in this stage as well. 
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5.3 Plateau, stage C— D 

This corresponds to the time interval t G [tctn] C [0,t*]. Again, we keep tc in some 
formulae as a symbolic reminder of the beginning of this asymptotic stage, 1 <^ tc* <^ ^2~^, 
and we define t^, = t^, precisely. The voltage E and the n-gate change on a time scale of 
the order — tc = 0{e2~^) along the upper systolic branch of the super-slow manifold 
until E{t^) = -E* and n{t^:) ~ N^,. The /i-gate remains close to its quasi-stationary value 
at fti 0. 

(a) Asymptotic solution. The general asymptotic results describing this stage are given 
by (30). Since the functions of the caricature model in the time interval t G [0,t*] are 
G{E) = k^i^Es — E), noo{E) = 1, g2{E) = fn{E) = En, the quadrature (30) evaluates 
to 

n = 1 - exp (-F„(t2 - t2,c)) , (49a) 
E = Es + f{l-exp (-F„(t2 - t2,c)))\ (49b) 
where ^2 c = ^^tc and lim ^2 c = 0- 

(b) Limit of the exact solution. The plateau stage occurs during the time interval 
t G [0,t*]. In order to compare the limit of the exact solution (35) to equations (49) we 
need to change to the super-slow time, ^2 = ^2t- Then the solution (35a) for n{t2) does not 
contain a small parameter and is readily comparable with the asymptotics, 

lim n{t2) = 1 - exp (-F„,(t2 - t2,c)) , (50) 

where we have taken into account that the initial value of n according to (44a) is n(t2,c) = 0- 

1 

For the voltage during this stage, we consider E from (35c) with a change to the 
super-slow time t2 = 62^, ^ 

E{e2-H2) =EBc{e2~h2) 

- - ( t - ( - 1 ) - it) 4'-^'' (') " ('^ - " - ■ 

(51) 

and look for the limit e — > 0, €2 — * at a fixed ^2- The terms denoted here by -E'bc(£2"^^2) 
are precisely those discussed in connection with the limit of the exact solution during stage 
B-C with the only difference that now we consider them in the super-slow time ^2- For a 
fixed t2, and small €2, expression (48) evaluates to 

lim EBc(e2~'t2) = E3. (52) 

£2^0 



•^Here and in similar cases further on, we imply that E and -Ebc are defined as functions of t, rather 
than t2 or any other stage-specific time variable. 
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Using expressions (45) derived in the previous stage the remaining u{-, ■)-function in (51) 
in the hmit e — becomes 

r ff. ,N , t2\ 1 -exp(fc3t2) exp(-(4-/)F„t2) 

hmu 4-/ e2F„-A;3, - = — . 53 

e^O \^ 62/ (4 - /) e2F„ - /C3 

The exponentially growing term is cancelled by exp(— A;3t2/e2) in (51) and therefore sub- 
stituting the above expressions in (51) and taking the limit 62 — > ultimately gives 

hm E{t2) = E3 + (1 - exp(-F„t2))' , (54) 

in full accord with the asymptotic result (49b) for the plateau stage of the AP, if we replace 
t2,c with its limit lim t2,c = 0- 

Before proceeding to the next stage, we comment on the position of the point D on the 
{n, E) plane, which is close to the end of the systolic branch of the super-slow manifold, 
the point (A^*, £'*). As stated earlier, we define point D by the exact condition E(t^) = E^, 
hence the condition n{t/) ^ is only approximate, since the AP trajectory follows the 
super-slow manifold only approximately, with precision 0(62). We shall see shortly that 
for the next stage it is important that n(t*) 7^ A^*. The sense of this inequality can be 
easily seen from (26a): we know that during the C-D stage E{t) decreases, G{E) > and 
g2{E) < 0, hence n > {—G{E)/g2{E))^^^ = Af{E) during the whole of that stage, which 
for t = t^, E = E^: gives n > M{E/) = N^. More accurately the value of n(t*) can be 
estimated using perturbation theory in €2 around the super-slow manifold n = Af{E); this 
would be a further distraction from our main goal, and we omit these formulae, like in the 
problem of determining and tj-. 



5.4 Repolarization, stage D— E 

During this stage the voltage E jumps on a time scale ~ 1 from the systolic to the diastolic 
branch of the super-slow manifold. The /i-gate changes swiftly on a time scale ~ e from 
its lower quasi-stationary value close to zero to its upper quasi-stationary value close to 
unity due to the discontinuity in the right-hand side of equation (33b). The slow n-gate 
remains approximately unchanged at n ~ A^*. The associated time interval is t G [tu, t^;] = 
[t^,tj] U [tf,t£;] where tf) = is the beginning of this stage, t| is the time of the inflexion 
point defined by E{t/) = E^ and Ie is the time of the end of the stage constrained by 
1 < -tf < £2"^- 

(a) Asymptotic solution. The asymptotics at this stage are given by quadrature (27). 
During this stage of the AP, the form of the caricature equations changes as the solution 
E(t) moves through the point E^. 

In the time interval t G [t*, tj] the relevant functions of the caricature model are G{E) = 
k2{E — E2) and g2{E) = g2 and therefore quadrature (27) evaluates to 

E=(^E2- ^^^) + (^E^ - {e2 - ^^^) j exp - t.)) , (55a) 



Asymptotic analysis of cardiac excitation (2007/04/04) 



30 



n = n{t/) = 1 - exp (-i^„(t2,* - h,c)) = const, (55b) 

with initial conditions and E{t^) = E^,. The function E{t) given by expression 

(55a) monotonically decreases since n (tj > iV, = {k2{E^ - E2)/g^Y'^. In the second time 
interval, t > t-\, the relevant functions of the caricature model are G{E) = ki[Ei — E) and 
g2{E) = 0, and quadrature (27) gives 

Eit) = (Et - El) exp {-hit - tt)) + El, ^^g^ 
n(t) = n{t/) = const, 

with initial conditions n{t^,) and E{t^) = E^. 

Finally, we note that the dynamics of h gate during this stage have a peculiarity due to 
the discontinuity of the right-hand side of (33b). The finite constraint h = 9{E^ — E) which 
according to (22b) is supposed to be approximately observed outside the AP upstroke, 
cannot be observed when E crosses the level E^, as this would mean a discontinuity in 
h{t). In fact, the jump of h{t) from 1 to happens gradually, of course. We can see from 
(33b) that this jump takes time ~ e. This violation of (22b) does not, however, in any 
way affect the dynamics of E and n, as E^ = Ej^ < E^ = E^ and therefore the factor 
h9{E — E*) in (33a) remains identically zero throughout this stage. 

(b) Limit of the exact solution. In order to compare the asymptotic and the exact 
solutions in a given AP stage we need to align them in time. The beginning of the present 
stage is tz? = t* = t2,*/e2, so we set 

t = €2-%^, + t (57) 

and then consider the limit e, e2 — > at a fixed t. As discussed earlier, we assume here 
that the limit of ^2,* is known (and finite). With this the exact solution in the interval 
t G [t*, tj-] becomes 

n{e2~^t2^^ + i) = 1 - exp(-F„t2,* - e2E„t), (58a) 

E{e2'H2,, + i) = {E, - w{e2~H2,/)) exp {k2 i) + u;(e2"42,* + i). (58b) 

We have immediately that in the limit €2 — > 0, expression (58a) coincides with (55b), since 
limt2,c = 0. Then, we have from (36b) 

lim w(^-^+i]=E2-^l V(-l)' e-'^"*^- = ^2 - (1 - exp(-F„t2,.)) . (59) 

According to (35a), n{t/) = 1 — exp(— F„e2t*) = 1 — exp(— F„t2,*)- Hence we have that the 
limit of (58b) coincides with (55a) since t = t — t^,. Analogously, substituting (57) in the 
exact solution (35) for the interval t > t^^ and taking the limit €2 — we reproduce the 
asymptotic solution (56) identically. 

Finally, the limit of the exact solution (35b) for the /i-gate as e — > results in /i = 
for t G [t*,tf] and h=l for t G in accordance with the asymptotic theory. 
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5.5 Recovery, stage E— F 

The voltage E and the n-gate move on a time scale ~ 62"^ close to the diastolic branch 
of the super-slow manifold. This continues forever as the system approaches it stable 
equilibrium. The associated time interval is t G [tE, +00) with tE as defined above. Since 
tE ~ ^2~^, we define t2^E = ^2tE which has a finite limit, coinciding with that of €2^1 and 

(a) Asymptotic solution. The asymptotic solution in the the recovery stage is given by 
quadrature (30). With the approximations of section 4.2, the functions of the caricature 
model in (30) are noo{E) = 0, fn{E) = Fn and N'~^{n) = Ei and so the asymptotic 
solution is 

n = n{tE) exp (-F„(t2 - h^E)) , (60a) 
E = El, (60b) 
h = l. (60c) 



(b) Limit of the exact solution. At a fixed super-slow time ^2 = ^2t, the limit €2, e ^ 
of the exact solution (35) is 

lim n(t2/e2) = (exp(i^„t2,t) - l) exp(-F„t2), (61a) 

hm 1(^2/62) = |m (Et - El) exp {-ki{t2 - t2,t)/e2) + ^1 = Ei, (61b) 
lim h{t2/e2)= lim (1 - (1 + exp(F;,t2,t/(ee2))) exp(-Fft,t2/(ee2))^i) = 1. (61c) 

Finally, we notice that according to (61a), n(t^) = 1 — exp(— F„t2,t)- Using this fact 
equation (61a) can be shown to be identical to equation (60a), so there is full agreement 
between (60) and (61). 



We have demonstrated that at every stage of the AP the explicit asymptotic solution 
coincides with the appropriate limit of the exact analytical solution of the caricature model 
(33). We conclude that the asymptotic theory is validated. 



6 Discussion 

6.1 Summary of results 

We have started from the classical Noble model [31] of cardiac Purkinje fibres, arguing 
that it is the simplest ionic model based on cardiac electrophysiology. Using numerical 
observations, we have postulated a system of axioms which allowed us to propose a reason- 
able parametric embedding (9) of the Noble model. We have also noted that some of the 
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features of the Noble model are rather peculiar in comparison with other cardiac models. 
This have lead us to propose the Archetypal Model (19) with the "generic" structure of 
modern cardiac models, which allows a simple parametric embedding (20) and which, in 
addition, is still a very accurate approximation of the Noble model and whose asymptotic 
limit coincides with that of the Noble model. Finally, we have obtained analytical solu- 
tions in quadratures, given by formulae (23), (24), (27) and (30), corresponding to the 
asymptotic limits in the embedded small parameters in the Archetypal Model. 

In that sense, we have achieved a fully analytical description of an action potential in 
a detailed ionic model of cardiac excitability. 

The accurate reproduction of the properties of the authentic ionic model necessitated 
a number of mathematical features of the parametric embedding used which made the 
standard singular perturbation approaches based on Tikhonov theorem inapplicable: 

• a large factor in front of only some terms in the right-hand side of the same equation; 

• non-analytical, perhaps even discontinuous, asymptotic limit of some right-hand 
sides, even though the original system is analytical, 

• non-isolated equilibria in the fast subsystem; 

• dynamic variables which change their character from fast to slow within one solution 
(remember in Tikhonov's theory, the roles of "fast" and "slow" variables are fixed); 

• the slow set may not even be a manifold, but may consist of pieces of different 
dimensionality (see Appendix B). 

All these features are related to each other and originate from the biophysics of excitable 
membranes, namely the fact that ionic gates work as nearly-perfect switches. Note that in 
some previous two-component simplified cardiac models, e.g. [33, 42? ], there are segments 
where null-clines of both variables are very close to each other. Perhaps, in an appropriate 
asymptotic embedding, these segments would be near continua of non-isolated equilibria 
in the fast subsystem, which is one of the key features mentioned above and which might 
be related to the success of those models. 

The asymptotic analysis we have done reproduces, in the limit e — 0, all five qual- 
itative phenomenological features of cardiac excitability, listed in section 1.1, which are 
inconsistent with FitzHugh-Nagumo type systems. Namely, 

1. Slow repolarization. The only fast part of a typical AP is the upstroke A-B, which 
goes from the upper edge of the red cross-hatched rectangle on figure 6(a) towards 
dashed blue line along the vertical axis. The other parts B-C-D-E-F all go along 
the dashed blue hne there, which is a set of the equilibria of the fast subsystem. In 
other words, all stages B-C-D-E-F are described in the slow subsystem (26). As seen 
in figure 6(b), this includes relatively fast parts B-C and D-E as well as relatively 
slow ones C-D and E-F, but these have different speeds due to the secondary small 
parameter €2- From the viewpoint of the main small parameter e they are all slow, 
i.e. much slower than the upstroke. 
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2. Slow subthreshold response. This corresponds to a trajectory starting from an 
initial voltage Eq < Em- For Eq above (below) the line h = H{E), this corresponds to 
the leftward (rightward) going trajectory in the red cross-hatched region of figure 6(a). 
That is, the only fast process, if any, following a subthreshold initial perturbation, is 
the relaxation of the h gate. The fast dynamics of E are not engaged as channels 
remained closed. Thus, any dynamics of E after such initial perturbation occur only 
on the slow time scale. 

3. Fast accommodation. In voltage-clamped conditions, i.e. when E{t) is prescribed, 
the fast subsystem (21) gives that h{t) tends towards H{E) 6{Eh — E), i.e. towards 
the blue dashed line on figure 6(a), and of course for E staying above E^ long 
enough, we have eventually h = 0. In other words, if E raises too slowly, the h 
gates have sufficient time to close, which prevents excitation. For this to happen, the 
(prescribed) dynamics of E should be slow compared to h dynamics, which are fast 
in terms of the e-embedding. 

4. Variable peak voltage. The peak voltage is the voltage at point B. From figure 6(a) 
it appears that this voltage is always very close to Ej^a as long as point A is above the 
line E = Em- However, according to the analysis of section 3.2.1, the peak voltage 
Eoo is determined via equation ho = J(-E'oo) — J{Eq), that is it does depend on initial 
condition. This paradox is explained in the end of section 3.2.1 and in Appendix B as 
a consequence of presence, in the Noble model and its descendant Archetypal Model, 
of a yet another small parameter a which is a factor of fh/gi- Namely, Ej^a ~ 
happens to be exponentially small in a. A conclusion follows from there that if the 
parameters in the Archetypal Model are changed in such a way that fh/di is not so 
small — e.g. by decreasing the value of Gj^ai the peak voltage will demonstrate a 
noticeable dependence on and Eq. This in fact happens in other cardiac models 
that are not so stiff as Noble model, i.e. in Courtemanche et al. [14] model of human 
atrium: as we have shown in [11], the phase portrait of its fast subsystem is very 
similar to figure 6(a), but less stiff and the peak voltage does indeed vary widely in 
the whole of {Em, E^a) range. 

5. Front dissipation. This feature requires the analysis of the spatially distributed 
version of the models and is therefore beyond the scope of this paper. However, the 
phenomenon of front dissipation has been a specific target of a number of our previous 
papers [6, 8, 11, 37]. In [6, 8] we considered a piecewise linear "caricature" version of 
the fast subsystem which allowed an explanation of front dissipation via establishing 
existence of a lower limit for the propagation speed, as opposed to FitzHugh-Nagumo 
type system which do not have such limit. In [11, 37] we demonstrated that the 
spatially extended version of fast system in [14], which as noted above is essentially 
identical to that of the AM, demonstrated the lower limit of the conduction velocity 
similar to that in the piecewiselinear caricature. Moreover, we have demonstrated 
that the understanding of conduction blocks based on this lower limit has a predictive 
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ability for wavebreaks in complicated spatiotemporal regimes in the Courtemanche 
et al. model. 

As no rigorous theory of slow-fast systems of non-Tikhonov type exists at present, 
we have formulated a "caricature" -style simplification of the AM, which is a less accurate 
approximation of the Noble model but has the same asymptotic structure as AM and admits 
exact solution. Using this exact solution we have been able to prove the asymptotic results 
in the particular case. 

6.2 Further directions 

We believe that the described procedure is generic and can be applied to typical detailed 
cardiac models including the complicated contemporary models. For instance, part of the 
same procedure has already been successfully applied to the human atrial tissue model 
of Courtemanche et al. [14] for which an analytical condition for propagation block in a 
re-entrant wave has been derived and a satisfactory quantitative agreement with results 
of direct numerical simulations have been demonstrated [6, 8, 37]. It is of even greater 
interest to investigate break-up and self-termination of AP fronts in models of ventricular 
myocytes such as [4, 23, 28, 39] since cases of ventricular fibrillation have more serious 
health consequences than those of atrial fibrillation. Another direction in which the pro- 
posed asymptotic description might be useful is the derivation of action potential restitution 
curves and conduction velocity dispersion curves for realistic cardiac models. These two 
curves are the most popular and widely available experimental characteristics of cardiac 
tissue upon which various interpretations of cardiac dynamics are based [36]. 

Finally, a simplified model, like the Archetypal Model or its caricature suggested here, 
can be a useful tool for large-scale numerics. Confidence in such a tool will increase if the 
simplified model has been derived by a controlled asymptotic procedure from a detailed 
model and preserves its predictive power. Such a model can also be useful for theoretical 
studies. For example, the difference between "slow over-threshold response" and "normal 
fast upstroke", discussed in Appendix C only appears in the asymptotic limit e and 
cannot be mathematically identified in the model at e = 1, even if the exact analytical 
solution like (35) is available. This difference may be important physiologically. E.g. it 
creates a possibility for "slow" and "fast" propagating waves in the system, a feature that 
has been observed in other models and in electro-physiological experiments as "Na" and 
"Ca" excitation waves [34, 35]. We believe that asymptotic analysis is the most adequate 
tool for mathematical description of this sort of "qualitative" phenomena, and cannot be 
replaced with numerical or even exact analytical solutions. 
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Figure 10: (Color online) (a) The function Q{E) defined by (62) in comparisson with rv?^(E) and 
hooiE). (b) Solutions {E,h,n) of the Noble model (4) (thick lines) and of the more accurate version of 
the Archetypal Model (66) (corresponding thin solid lines). 



Appendices 

A A more accurate archetypal model 

As one can see from figure 4(a), the AP of the AM (19) is a bit longer, and its repolarization 
is a bit slower than that of the Noble model (4). Although this difference may seem 
relatively minor, it is in fact surprisingly large, considering that the small parameter e 
used to derive equations (19) from equations (4) is related to small quantities in (4) of the 
order of 10^^. Thus, we would expect an accuracy of the order of 1% in all results, which 
is clearly not the case in figure 4(a). The reason for this is that the asymptotic structure 
of Noble model (4) is even more complicated than that summarized in the Axioms I- VII. 
Here we argue that the observed discrepancy is mainly due to a deficiency of Axiom VII. 
However, we believe that the complication in question is idiosyncratic for Noble model and 
is not actually observed in later more realistic models. Still, in order to demonstrate the 
validity of our approach, we show here how the AM can be improved by an appropriate 
correction. 
f:0070 

Figure 2(b) shows that the function S{E) is relatively small. However, it appears in 
(15) multiplied by the large factor gi{E) / fh{E) and thus the values of the product 

Q{E)^g,{E)S{E)/ME), (62) 

are in fact of the order unity over a significant range of voltages as shown in figure 10(a). 
It is then clear that neglecting this essentially non-zero term in equation (15) leads to 
the low accuracy of the AM (19). At first glance, this means that in the Noble model 
(4), during the repolarization phase of the AP, gate h is, in fact, not fast compared to 
other two variables. This would mean that reduction to a two-variable model in that 
region is not possible. However, as we have discussed in Introduction after the definition 
of embedding, a replacement of 1 with a small parameter can, in fact, be a reasonable 
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embedding, and its quality can be assessed by a comparison of the numerical solutions of 
the original and the embedded problem. Such comparison is provided in figure 3 and shows 
that the embedding (19) is indeed reasonable although of a not very high accuracy. Thus 
we suppose that a higher accuracy can be achieved by taking a first-order approximation of 
the term Q in the parameter e, rather than zero-order as in the other terms. Let us denote 
the small parameter associated with the function Q{E) hj fi, fi > 0. Then the error term 
after the first-order approximation in fi will be 0{^'^). The error term after the zero-order 
approximation in e will be, of course 0(e). If we want these two kinds of error terms to be 
of the same order, we must therefore consider /i = e^^^. These arguments are formalized 
by the following improved version of Axiom VII, 

Axiom Vila. S{E; e) = e^/^S{E) + 0(e), where S{E) ^ S{E). 

Besides, for technical reasons the following, stronger version of Axiom VI will be more 
convenient for us: 

Axiom Via. m?^{E] e) h^{E] e) = eW{E) + O(e^), for the same W{E) as in Axiom VI. 
Substituting Axioms Via and Vila into (15), we get 

:E.„„,,_,.|(|„,)^^,o„. (63, 

From here we deduce that dE/dt = giW + 0(e^/^). Substituting this into the right-hand 
side of (63) we obtain 

^ = amE) - eV^llll S{E)g,W + 0(e). (64) 

Thus, after discarding 0(e) and putting e = 1, we arrive at the following variant of (16), 

^ = (^?i(E) W{E) + g^iE) + g,{E)) (l - Q{E)) , (65a) 
h = H{E)e{Eh- E), (65b) 

— = fn{E) {n^{E)-n), (65c) 
dt 

where the first and third equations are satisfied with accuracy 0(e) and the second equation 
is satisfied only with accuracy 0(e^/^). Here Q{E) = S{E) gi{E) / fh{E) which according 
to Axiom Vila should be close to Q{E). 

We see that in terms of the slow subsystem, the modification simply amounts to mul- 
tiplying the right-hand side of equation (65a) by a known function of E. By analogy, it is 
straightforward to propose an improved version of the AM (19), 

^ = {g^{E) M{E) 9{E ~ E^) h + g,{E) W{E) + g,{E) + g^{E)) (1 - Q{E)) , 

^ = UE) {H{E) 9{E, -E)-h), (66) 

-^ = fniE) ME)-n). 
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This improved AM gives solutions very close to those of the Noble system (4), as illustrated 
in figure 10(b). 

Notice that the asymptotic structure of the improved AM (66) is exactly the same as 
that of (19). Indeed, the difference amounts to redefining the functions gi{E), g2{E), QsiE) 
via a factor 1 — Q, thus the asymptotic theory discussed in section 3.2 is equally applicable 
to both systems. 



B The high-excitability embedding 

For the numerical values of the parameters corresponding to healthy tissue, the voltage 
upstroke at the beginning of the AP is, in fact, a much faster variable than the /i-gate. 
Indeed, the voltage speed constant gNa/Cu = 33| is large compared to the typical values 
of the /i-gate speed function, max fh{E) ~ 1. This speed difference is unaccounted for 

by Axioms I-VII(a) where these two quantities have the same asymptotic order. How- 
ever, some features of a typical AP solution depend on the ratio of these quantities in an 
exponential way. 

To take this feature into account, here we construct an improved embedding which 
takes the E-vs.-h speed difference into account. We use an additional embedding with one 
more artificial small parameter cr > 0. Formally, we replace Axiom I with a new Axiom, 

Axiom la. gNa{e, a) = a'^e^'^gNa- 

This, of course, does not affect the slow-time subsystem (16). The fast-time subsystem 
(10) now depends on the new parameter cr and becomes 



^ = a-'g^{E)M{E)9{E-Ejh, 
^ = ME) {H{E)e{E,-E)-h) 



(67) 



where the trivial equation for n is omitted as the n-gate neither changes nor matters on 
time scales of interest here. The equations (67) appear to be a standard fast-slow Tikhonov 
system. However, this is deceptive since the specific properties of the right-hand sides lead 
to a number of nonstandard features. Let us consider the limit a — > in this system. In 
the super-fast time s = T/a, the /i-gate is a first integral, dh/ds = 0, and E satisfies the 
super-fast equation, 

dE 

— = gi{E) M{E) e{E - E„) h, (68) 
ds 

depending on /i as a parameter. The equilibria of this system for which the the right-hand 
side vanishes, form an unusual set consisting of the lines /i = and E = Ejsfa and the semi- 
stripe {{E,h)} = {—oo, Em) X [0,1]. Hence, the slow set is not a manifold, but consists 
of pieces of different dimensionalities. This feature is even "more non-Tikhonov" than the 
many non-standard properties observed in the main embedding. One consequence is that 
the slow-time subsystem of (67) changes form as the trajectory passes through the various 
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pieces of the slow set. This dependence is substantiah even the dimensionahty of the slow 
system changes. For instance, for a typical AP solution in the beginning of the plateau, 
the trajectory crosses the E = E^a piece, which can be parameterized by h. Then we have 
a one-dimensional slow subsystem of (68) (which is in the "fast time" T) in the form of an 
equation for /i, 

^ = -UE^a) h. (69) 

Then, during the later part of the plateau, which proceeds along the one-dimensional piece 
/i = 0, the slow subsystem becomes "zero- dimensional" in the sense that all right-hand 
sides vanish and all trajectories are fixed points. This corresponds to the fact that the 
movement along this piece occurs slowly in terms of the parameter e, i.e. is infinitely slow 
not only in terms of time s but in terms of time T too. Finally, during the repolarization 
phase of the AP when the voltage E drops below Em, the slow subsystem is on the two- 
dimensional piece but is in fact foliated to one- dimensional pieces, as the right-hand side 
of the equation for E vanishes and the voltage ii^ is a first integral, 

dE 

^ = h{E) {H{E)e{E,-E)-h). 

Since the equations of this "high-excitability" embedding are at most one-dimensional 
systems, obviously all of them can be solved in quadratures. This embedding is rather 
instructive since it demonstrates that a non-Tikhonov embedding might lead to rather 
untypical consequences. It might also be useful for a number of applications, especially 
when healthy, well-excitable tissues are concerned. However, the most important applica- 
tions are those related to the failure of excitation, or of excitation propagation in the case 
of spatially-extended systems. These processes are observed, however, exactly when the 
excitability, represented here by formal parameter a~^, is not so high. 



C Asymptotic synthesis 

We assume, for simplicity, that the initial values of the gating variables are given by /iq = 1, 
TT-o = 0. Then, depending on the initial trans-membrane voltage Eq, there exist three types 
of solutions of (19) as visualized by the phase portraits in figure 6 and described below: 

C.l Sub-threshold response 

If the initial value of E is less that the threshold value of the super-slow subsystem (28) 
i.e. Eq < E2, the voltage decays towards its global equilibrium Ei. 

Relaxation of the h-gate. The /i-gate relaxes on a time scale ~ e towards its quasi- stationary 
value h ~ H{E) 9{Eh — E) according to (23). The variables E and n remain close to 
their original values of E ^ Eq, n ^ uq = 0. 
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Relaxation of the voltage E. The voltage decays on a time scale ~ 1 towards its equilibrium 
value El according to (27) (with n = 0). The /i-gate remains close to its quasi-stationary 
value except, possibly, swiftly moving on a time scale ~ e along its discontinuity as E 
passes through Eh- The slow ra-gate remains approximately at n no = 0. 

C.2 Slow over-threshold response 

If the initial value of the voltage is bigger than the threshold value of the super-slow 
subsystem (28) but smaller than E^ i.e. E2 < Eq < E^ then the slow subsystem alone is 
sufficient to describe the AP evolution and the fast Na current is not involved. The voltage 
makes a relatively small excursion towards the upper systolic branch of the super-slow 
manifold and approaches it from below: 

Relaxation of the h-gate. The h-gate relaxes on a time scale ~ e in the same way as in the 
case of sub-threshold response, except this time the quasi-stationary value of h is zero 
since E2 > Eh. 

Rise of the voltage E. The voltage increases on time scale ~ 1 towards E3 i.e. towards the 
upper part of the systolic branch of the super-slow manifold according to (27) (with 
n = 0). The slow ra-gate remains approximately unchanged at n ~ no = 0. 

Plateau. The variables E and n move on a time scale ~ e2~^ along the upper systolic 
branch of the super-slow manifold n = //{E), according to (28) until they reach the 
point [E^:, N^,). The /i-gate remains close to its quasi-stationary value of ~ 0. 

Repolarization. The voltage E jumps on a time scale ~ 1 according to (27), towards 
the diastolic branch of the super-slow manifold. The /i-gate remains close to its quasi- 
stationary value except, possibly, for a swift movement along the discontinuity of H{E)6{Eh~ 
E) aX E = Eh on a time scale ~ e. The slow n-gate remains approximately at n ~ A^*. 

Recovery. The variables E and n move on a time scale ~ e2~^ along the diastolic branch 
of the super-slow manifold n = JV{E), according to (28). This continues forever, with 
{E,n) asymptotically approaching the true equilibrium of the system {Ei,0). Gate h 
stays close to its quasi-stationary value h ^ 0. 

C.3 Normal fast-upstroke action potential 

If the initial value of the voltage exceeds the threshold of the primary fast-time system (21) 
i.e. Eo > Em the fast Na current is activated and a normal fast-upstroke AP, is initiated: 
Fast upstroke. The variables h and E change together on a time scale ~ e according to 
(24) from the point {E, h) = {Eq, 1) asymptotically (in T — > 00) approaching the point 
{E,h) = {Eoo{Eq,1),0). The slow n-gate remains approximately unchanged at n ^ 
no = 0. 

Post- overshoot drop of the voltage E. The voltage E descends, inasmuch as Eoo{Eq,1) > 
E-^, on a time scale ~ 1 towards its higher equilibrium value E^ according to (27) (with 
n = 0). Variable h remains close to its quasi-stationary value of /i ~ 0. The slow n-gate 
remains approximately unchanged at n ~ no = 0. 



Asymptotic analysis of cardiac excitation (2007/04/04) 



40 



Plateau, repolarization and recovery stages follow which are similar to the corresponding 
stages in the case of slow over-threshold response. 
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